************************************************************************************************
* This replication file maps the profitability at the municipal level after a price reduction. *
************************************************************************************************
* 1. Run the file which re-calibrates the parameters to correspond to the price reduction selected in 07_01_price_reduction_loop.do.
* 2. Calculate output and profits and aggregate it to the municipal level.
* 3. Map the results.

*************************
* 0 Set path to files *
*************************
clear all
capture log close
set more off
set trace off 
mata: mata set matastrict off
cd $main_directory


**************************************************************************************************************************************
* 1. Run the file which re-calibrates the parameters to correspond to the price reduction selected in 07_01_price_reduction_loop.do. *
**************************************************************************************************************************************

//global reduction_level_Q1 0 //This sets the % price reduction for real estate transactions if you would like to set it manually.
//global reduction_level_Q2 0 //This sets the % price reduction for other transactions if you would like to set it manually.
do 07_02_price_reduction_function.do //This runs the code which re-calibrates the parameters to reflect the price reduction

****************************************************************************
* 2. Calculate output and profits and aggregate it to the municipal level. *
****************************************************************************

//Predict output
sort market_ID
mata:
 X    = st_data(., "$controls_total_Q1") //This selects the corresponding control variables.
    nx    = rows(X)
    X    = X,J(nx,1,1)
	Results=spdemand_log(X,betas_Q1_new)
	n_firms=rows(Results)
    st_numscalar("n_firms",n_firms)
	firm_ID    = st_data(., "firm_ID") 
	Prediction=J(nx,1,.)
	for(i=1; i<=nx; i++) { //i goes over the market-notary combinations
    for(j=1; j<=n_firms; j++) {
	if (firm_ID[i]==Results[j,1]) {
		Prediction[i]=Results[j,2]
	}
	}
	}
prediction_exp=exp(Prediction)
st_addvar("double","Q1_IV")
st_store(.,"Q1_IV",prediction_exp)
end

sort market_ID
mata:
 X    = st_data(., "$controls_total_Q2") //This selects the corresponding control variables.
    nx    = rows(X)
    X    = X,J(nx,1,1)
	Results=spdemand_log(X,betas_Q2_new)
	n_firms=rows(Results)
    st_numscalar("n_firms",n_firms)
	firm_ID    = st_data(., "firm_ID") 
	Prediction=J(nx,1,.)
	for(i=1; i<=nx; i++) { //i goes over the market-notary combinations
    for(j=1; j<=n_firms; j++) {
	if (firm_ID[i]==Results[j,1]) {
		Prediction[i]=Results[j,2]
	}
	}
	}
prediction_exp=exp(Prediction)
st_addvar("double","Q2_IV")
st_store(.,"Q2_IV",prediction_exp)
end

//We then calculate costs and profits (this requires only 1 observation per notary).
keep if common==1
sort firm_ID
gen L2_Q=exp(-aL)*(alpha*Q1_IV+(1-alpha)*Q2_IV)^(1/nu)
gen M2_Q=exp(-aM/phi)*(alpha*Q1_IV+(1-alpha)*Q2_IV)^(1/(nu*phi))
gen VC_Q=L2_Q+M2_Q
gen v_profit_individual_n=p_Q1_new*Q1_IV+p_Q2_new*Q2_IV-VC_Q
gen profit_individual=v_profit_individual_n-fixed_cost*NotPerOff
count if profit_individual<0
keep nis_not profit_individual
bysort nis_not: egen mean_mun_profit=mean(profit_individual)
bysort nis_not mean_mun_profit: gen dup = cond(_N==1,0,_n)
drop if dup>1
drop profit_individual
//Store the information
save "Estimates\Graphs\profitability_municipality_`=$reduction_level_Q1'_`=$reduction_level_Q2'.dta", replace

***********************
* 3. Map the results. *
***********************

clear
//The files used in the analysis were generated from the original layers of spatial data using the commented out commands:
//ssc install spmap
//ssc instal shp2dta
//Turn the shapefile (GIS) file into a stata file
//shp2dta using CENSUS_UNIT_RG_01M_2011, database(units_db) coordinates(units_coord) genid(id)
//Open the file with the geographic data
use "Data_Belgium\GIS_data_census_units\units_coord.dta"
geo2xy _Y _X, replace proj(mercator)
save "Data_Belgium\GIS_data_census_units\geo2xy_coor_proj.dta", replace
use "Data_Belgium\GIS_data_census_units\units_db.dta", clear
//Generate a country identifier
gen country=substr(CENSUS_ID,1,2)
//Keep only the observations in Belgium
keep if country=="BE"
//Generate the NIS variable
gen nis_not=substr(CENSUS_ID,7,5)
destring nis_not, replace
//Merge with data on relevant notary-level variable
merge 1:1 nis_not using "Estimates\Graphs\profitability_municipality_`=$reduction_level_Q1'_`=$reduction_level_Q2'.dta"
drop _merge
//Install color palette aid
//ssc install palettes 
//ssc install colrspace
colorpalette greys, n(5) nograph
local colors `r(p)'
if $reduction_level_Q1 ==0 & $reduction_level_Q2 == 0 {
spmap mean_mun_profit using Data_Belgium\GIS_data_census_units\geo2xy_coor_proj, id(id) fcolor("`colors'") clmethod(custom) clbreaks(-100 0 100 200 300 1000) legend(size(large)) legend(label(1 "no entrant"))
}
else if $reduction_level_Q1 !=0 & $reduction_level_Q2 !=0{
spmap mean_mun_profit using Data_Belgium\GIS_data_census_units\geo2xy_coor_proj, id(id) fcolor("`colors'") clmethod(custom) clbreaks(-100 0 100 200 300 1000) legend(off)
}
graph save Estimates\Graphs\profitability_municipality_`=$reduction_level_Q1'_`=$reduction_level_Q2', replace
//Paper map
graph export "Estimates\Graphs\profitability_municipality_`=$reduction_level_Q1'_`=$reduction_level_Q2'.eps", replace
